#pragma once
#include <algorithm>
#include "gtest/gtest.h"
#include "DataStructure/Common/Id.h"
#include "DataStructure/Input/Solver/Subcase.h"
#include "DataStructure/Input/Property/PSOLIDData.h"
#include "DBManager/DBServiceFactory.h"
#include "DBManager/IService/Input/Coordinate/ICORD2RService.h"
#include "DBManager/IService/Input/Material/IIsotropicMaterialService.h"
#include "DBManager/IService/Input/Material/IAnisotropic3DMaterialService.h"
#include "SGFem/Pre/Assembly/Assembly.h"
#include "SGFem/Pre/Assembly/CoordinateProcessor.h"
#include "../../Factory/ElementCalculatorFactory.h"
#include "../SolidMaterialParser.h"
#include "../Tetra4ElementCalculator.h"

using SG::Algebra::Matrixd;
using SG::Algebra::Point;
using SG::Algebra::PointArray;
using SG::DataStructure::Common::Id;
using SG::DataStructure::Common::IdArray;
using SG::DataStructure::Common::Real;
using namespace SG::DataStructure::FEM;
using namespace SG::DBManager;
using SG::FEMSolver::AssemblyEigen;
using SG::FEMSolver::CoordinateProcessor;
using SG::FEMSolver::Element::ElementCalculatorFactory;
using SG::FEMSolver::Element::SolidMaterialParser;
using SG::FEMSolver::Element::Tetra4ElementCalculator;

// ElementCalculatorFactory 是否正确创建 Tetra4 单元计算类
TEST (ElementCalculatorFactoryTest, CreateTetra4ElementCalculator)
{
    ElementCalculatorFactory::Register (ElementTypes::ELETYPE_Tetra4, Tetra4ElementCalculator::Create);

    // Test data
    const auto calculatorPtr = ElementCalculatorFactory::Create (ElementTypes::ELETYPE_Tetra4);
    const auto tetra4CalculatorPtr = dynamic_cast<Tetra4ElementCalculator*> (calculatorPtr.get ());

    // Check
    ASSERT_TRUE (nullptr != tetra4CalculatorPtr);
}

extern std::shared_ptr<DBServiceFactory> DB_SERVICEFACTORY;
class Tetra4CalculatorTest : public ::testing::Test
{
protected:
    static void SetUpTestSuite ()
    {
        DB_SERVICEFACTORY->ClearDatabase ();

        // 坐标系
        const auto& cord2rService = DB_SERVICEFACTORY->get<ICORD2RService> ();
        CORD2RData rectangularCoordSys{100, 0, Point{1.0, 2.0, 3.0}, Point{2.0, 4.0, 1.0}, Point{5.0, 1.0, 1.0}};
        cord2rService->Add (rectangularCoordSys);
        rectangularCoordSys = {101, 100, Point{4.0, 3.0, 7.0}, Point{2.0, 6.0, 1.0}, Point{5.0, 1.0, 4.0}};
        cord2rService->Add (rectangularCoordSys);

        // 材料
        IsotropicMaterial material1;
        material1.m_id = 1;
        material1.m_E  = 210000.0;
        material1.m_nu = 0.3;
        material1.m_G  = 0.5 * material1.m_E / (1.0 + material1.m_nu);
        material1.m_rho = 7.85E-9;
        DB_SERVICEFACTORY->get<IIsotropicMaterialService> ()->Add (material1);

        Anisotropic3DMaterial material9;
        material9.m_id    = 9;
        material9.m_G[0]  = 1.067E7;
        material9.m_G[1]  = 6910000.0;
        material9.m_G[2]  = 6890000.0;
        material9.m_G[3]  = 1250000.0;
        material9.m_G[4]  = -90000.0;
        material9.m_G[5]  = 860000.0;
        material9.m_G[6]  = 1.903E7;
        material9.m_G[7]  = 6900000.0;
        material9.m_G[8]  = 1270000.0;
        material9.m_G[9]  = -120000.0;
        material9.m_G[10] = 830000.0;
        material9.m_G[11] = 1.428E7;
        material9.m_G[12] = 1230000.0;
        material9.m_G[13] = -150000.0;
        material9.m_G[14] = 820000.0;
        material9.m_G[15] = 5270000.0;
        material9.m_G[16] = 210000.0;
        material9.m_G[17] = -110000.0;
        material9.m_G[18] = 4630000.0;
        material9.m_G[19] = 520000.0;
        material9.m_G[20] = 5720000.0;
        DB_SERVICEFACTORY->get<IAnisotropic3DMaterialService> ()->Add (material9);

        // 生成节点排序
        auto& assembly = AssemblyEigen::GetInstance ();
        assembly.Clear ();
        Subcase subcase;
        assembly.Initialize (DB_SERVICEFACTORY, subcase, {1});

        // 计算坐标系原点和基矢量
        auto& coordProcessor = CoordinateProcessor::GetInstance ();
        coordProcessor.Clear ();
        coordProcessor.Initialize (DB_SERVICEFACTORY);
    }

    void SetUp () override
    {
        dispCoordIdSet = {0, 0, 0, 0};
        pointSet = {Point{0.0, 0.0, 0.0}, Point{0.9, 0.1, 0.0}, Point{0.2, 0.8, 0.1}, Point{0.3, 0.2, 0.7}};

        rho = 7.85E-9;
        D = Matrixd (6, 6);
        property.m_mId = 1;

        Tetra4ElementCalculator::CoordTrans (pointSet, _OUT T, _OUT localCoord);
        TEle = Matrixd (12, 12);
        for (std::size_t iNode{0}; iNode < 4; ++iNode)
        {
            const std::size_t i{iNode * 3};
            TEle.BlockFill (i, i, T);
        }

        eleDisplacement = Matrixd (12, 1, {
        0.0,                   0.0,                    0.0,
        0.0,                   0.0,                    0.0,
        0.0,                  -3.339067416526824E-6,   0.0,
        6.892407445945625E-5,  1.5768532222622028E-4,  7.64509518981474E-5
        });
        eleDisplacement = TEle.Transpose () * eleDisplacement;

        eleCalculator.Initialize (ElementTypes::ELETYPE_Tetra4);

        // 小量放大系数
        factor = 1.0E12;
    }

    IdArray dispCoordIdSet;
    PointArray pointSet;
    Real rho;
    Matrixd D;
    PSOLIDData property;
    Matrixd T;
    Matrixd localCoord;
    Matrixd TEle;
    Matrixd eleDisplacement;
    Tetra4ElementCalculator eleCalculator;
    Real factor;
};

TEST_F (Tetra4CalculatorTest, CoordTrans)
{
    // Test data
    Matrixd TTest (3, 3);
    Matrixd localCoordTest (3, 4);
    Tetra4ElementCalculator::CoordTrans (pointSet, _OUT TTest, _OUT localCoordTest);

    // Target data
    const Matrixd TTarget (3, 3, {
    -1.97055854E-01,  7.54723246E-01,  6.25749001E-01,
     7.09673161E-01, -3.30555379E-01,  6.22171316E-01,
     6.76411853E-01,  5.66679772E-01, -4.70468959E-01
    });
    const Matrixd localCoordTarget (3, 4, {
    0.00000000E+00, -1.06382953E-01,  5.95968543E-01,  5.56306173E-01,
    0.00000000E+00,  6.46195383E-01, -5.68316771E-02,  5.56981738E-01,
    0.00000000E+00,  6.25391233E-01,  5.75839957E-01, -1.71693074E-02
    });

    // Check
    ASSERT_TRUE (TTarget == TTest);
    ASSERT_TRUE (localCoordTarget == localCoordTest);
}

TEST_F (Tetra4CalculatorTest, ComputeStiffness)
{
    // Test data
    Matrixd TMat (3, 3);
    SolidMaterialParser::ComputeElasticMatrix (DB_SERVICEFACTORY, property.m_mId, TMat, _OUT D, _OUT rho);

    Matrixd kgTest (12, 12);
    eleCalculator.ComputeStiffness (localCoord, D, _OUT kgTest);
    kgTest = TEle * kgTest * TEle.Transpose ();

    // Target data
    const Matrixd kgTarget (12, 12, {
     3.15651822E+04,  1.46234818E+04,  1.19028340E+04, -2.23858300E+04, -4.33603239E+03, -1.27530364E+03, -2.85668016E+03, -1.20020243E+04,  3.75506073E+03, -6.32267206E+03,  1.71457490E+03, -1.43825911E+04,
     1.46234818E+04,  2.83414980E+04,  1.06629555E+04, -8.37449393E+03, -6.70242915E+02,  4.74696356E+03, -7.29048583E+03, -2.44305668E+04, -3.50708502E+03,  1.04149798E+03, -3.24068826E+03, -1.19028340E+04,
     1.19028340E+04,  1.06629555E+04,  2.39204453E+04, -5.31376518E+03,  4.07388664E+03,  9.38056680E+02,  3.08198381E+03, -7.54554656E+03, -1.51761134E+03, -9.67105263E+03, -7.19129555E+03, -2.33408907E+04,
    -2.23858300E+04, -8.37449393E+03, -5.31376518E+03,  3.04004049E+04, -4.20850202E+03, -7.65182186E+03, -4.86315789E+03,  1.46801619E+04, -3.04655870E+03, -3.15141700E+03, -2.09716599E+03,  1.60121457E+04,
    -4.33603239E+03, -6.70242915E+02,  4.07388664E+03, -4.20850202E+03,  1.05977733E+04,  1.55870445E+03,  9.96862348E+03, -7.09493927E+03, -2.86943320E+03, -1.42408907E+03, -2.83259109E+03, -2.76315789E+03,
    -1.27530364E+03,  4.74696356E+03,  9.38056680E+02, -7.65182186E+03,  1.55870445E+03,  1.25744939E+04, -1.70040486E+03, -4.88866397E+03, -5.95141700E+01,  1.06275304E+04, -1.41700405E+03, -1.34530364E+04,
    -2.85668016E+03, -7.29048583E+03,  3.08198381E+03, -4.86315789E+03,  9.96862348E+03, -1.70040486E+03,  1.23718623E+04, -3.12449393E+03,  7.43927126E+02, -4.65202429E+03,  4.46356275E+02, -2.12550607E+03,
    -1.20020243E+04, -2.44305668E+04, -7.54554656E+03,  1.46801619E+04, -7.09493927E+03, -4.88866397E+03, -3.12449393E+03,  4.01451417E+04, -6.69534413E+03,  4.46356275E+02, -8.61963563E+03,  1.91295547E+04,
     3.75506073E+03, -3.50708502E+03, -1.51761134E+03, -3.04655870E+03, -2.86943320E+03, -5.95141700E+01,  7.43927126E+02, -6.69534413E+03,  1.36188259E+04, -1.45242915E+03,  1.30718623E+04, -1.20417004E+04,
    -6.32267206E+03,  1.04149798E+03, -9.67105263E+03, -3.15141700E+03, -1.42408907E+03,  1.06275304E+04, -4.65202429E+03,  4.46356275E+02, -1.45242915E+03,  1.41261134E+04, -6.37651822E+01,  4.95951417E+02,
     1.71457490E+03, -3.24068826E+03, -7.19129555E+03, -2.09716599E+03, -2.83259109E+03, -1.41700405E+03,  4.46356275E+02, -8.61963563E+03,  1.30718623E+04, -6.37651822E+01,  1.46929150E+04, -4.46356275E+03,
    -1.43825911E+04, -1.19028340E+04, -2.33408907E+04,  1.60121457E+04, -2.76315789E+03, -1.34530364E+04, -2.12550607E+03,  1.91295547E+04, -1.20417004E+04,  4.95951417E+02, -4.46356275E+03,  4.88356275E+04
    });

    // Check
    ASSERT_TRUE (kgTarget == kgTest);
}

TEST_F (Tetra4CalculatorTest, ComputeLumpedMass)
{
    // Test data
    Matrixd lumpedMTest (12, 12);
    Tetra4ElementCalculator::ComputeLumpedMass (localCoord, rho, _OUT lumpedMTest);

    // Target data
    Matrixd lumpedMTarget (12, 12);
    for (std::size_t i{0}; i < 12; ++i)
    {
        lumpedMTarget (i, i) = 1.55364583E-10;
    }

    // Check
    ASSERT_TRUE (lumpedMTarget * factor == lumpedMTest * factor);
}

TEST_F (Tetra4CalculatorTest, ComputeCoupledMass)
{
    // Test data
    auto coupledMTest = eleCalculator.ComputeCoupledMass (localCoord, rho);
    if (std::any_of (dispCoordIdSet.begin (), dispCoordIdSet.end (), [] (Id id) {return id > 0;}))
    {
        coupledMTest = TEle * coupledMTest * TEle.Transpose ();
    }

    // Target data
    const Matrixd coupledMTarget (12, 12, {
     6.21458333E-11,  9.69352280E-27, -3.23117427E-27,  3.10729167E-11,  0.00000000E+00, -1.61558713E-27,  3.10729167E-11,  4.84676140E-27, -1.61558713E-27,  3.10729167E-11,  0.00000000E+00, -1.61558713E-27,
     9.69352280E-27,  6.21458333E-11, -9.69352280E-27,  0.00000000E+00,  3.10729167E-11,  0.00000000E+00,  4.84676140E-27,  3.10729167E-11, -4.84676140E-27,  0.00000000E+00,  3.10729167E-11,  0.00000000E+00,
    -3.23117427E-27, -9.69352280E-27,  6.21458333E-11,  0.00000000E+00, -1.61558713E-27,  3.10729167E-11,  0.00000000E+00, -4.84676140E-27,  3.10729167E-11,  0.00000000E+00, -1.61558713E-27,  3.10729167E-11,
     3.10729167E-11,  0.00000000E+00,  0.00000000E+00,  6.21458333E-11,  3.23117427E-27,  0.00000000E+00,  3.10729167E-11,  3.23117427E-27,  1.61558713E-27,  3.10729167E-11,  3.23117427E-27,  1.61558713E-27,
     0.00000000E+00,  3.10729167E-11, -1.61558713E-27,  3.23117427E-27,  6.21458333E-11, -3.23117427E-27,  3.23117427E-27,  3.10729167E-11, -4.84676140E-27,  3.23117427E-27,  3.10729167E-11, -4.84676140E-27,
    -1.61558713E-27,  0.00000000E+00,  3.10729167E-11,  0.00000000E+00, -3.23117427E-27,  6.21458333E-11,  1.61558713E-27, -1.61558713E-27,  3.10729167E-11,  1.61558713E-27, -1.61558713E-27,  3.10729167E-11,
     3.10729167E-11,  4.84676140E-27,  0.00000000E+00,  3.10729167E-11,  3.23117427E-27,  1.61558713E-27,  6.21458333E-11,  3.23117427E-27,  0.00000000E+00,  3.10729167E-11,  3.23117427E-27,  1.61558713E-27,
     4.84676140E-27,  3.10729167E-11, -4.84676140E-27,  3.23117427E-27,  3.10729167E-11, -1.61558713E-27,  3.23117427E-27,  6.21458333E-11, -3.23117427E-27,  3.23117427E-27,  3.10729167E-11, -4.84676140E-27,
    -1.61558713E-27, -4.84676140E-27,  3.10729167E-11,  1.61558713E-27, -4.84676140E-27,  3.10729167E-11,  0.00000000E+00, -3.23117427E-27,  6.21458333E-11,  1.61558713E-27, -1.61558713E-27,  3.10729167E-11,
     3.10729167E-11,  0.00000000E+00,  0.00000000E+00,  3.10729167E-11,  3.23117427E-27,  1.61558713E-27,  3.10729167E-11,  3.23117427E-27,  1.61558713E-27,  6.21458333E-11,  3.23117427E-27,  0.00000000E+00,
     0.00000000E+00,  3.10729167E-11, -1.61558713E-27,  3.23117427E-27,  3.10729167E-11, -1.61558713E-27,  3.23117427E-27,  3.10729167E-11, -1.61558713E-27,  3.23117427E-27,  6.21458333E-11, -3.23117427E-27,
    -1.61558713E-27,  0.00000000E+00,  3.10729167E-11,  1.61558713E-27, -4.84676140E-27,  3.10729167E-11,  1.61558713E-27, -4.84676140E-27,  3.10729167E-11,  0.00000000E+00, -3.23117427E-27,  6.21458333E-11
    });

    // Check
    ASSERT_TRUE (coupledMTarget * factor == coupledMTest * factor);
}

TEST_F (Tetra4CalculatorTest, ComputeCoupledMassWithLeftHandedGridOrder)
{
    // Input data
    pointSet = {Point{0.0, 0.0, 0.0}, Point{0.2, 0.8, 0.1}, Point{0.9, 0.1, 0.0}, Point{0.3, 0.2, 0.7}};

    // Test data
    Tetra4ElementCalculator::CoordTrans (pointSet, _OUT T, _OUT localCoord);
    auto coupledMTest = eleCalculator.ComputeCoupledMass (localCoord, rho);
    if (std::any_of (dispCoordIdSet.begin (), dispCoordIdSet.end (), [] (Id id) {return id > 0;}))
    {
        for (std::size_t iNode{0}; iNode < 4; ++iNode)
        {
            const std::size_t i{iNode * 3};
            TEle.BlockFill (i, i, T);
        }
        coupledMTest = TEle * coupledMTest * TEle.Transpose ();
    }

    // Target data
    const Matrixd coupledMTarget (12, 12, {
     6.21458333E-11,  2.26182199E-26, -1.61558713E-26,  3.10729167E-11,  1.13091099E-26, -8.07793567E-27,  3.10729167E-11,  1.13091099E-26, -8.07793567E-27,  3.10729167E-11,  1.13091099E-26, -8.07793567E-27,
     2.26182199E-26,  6.21458333E-11, -1.29246971E-26,  1.13091099E-26,  3.10729167E-11, -6.46234854E-27,  1.13091099E-26,  3.10729167E-11, -6.46234854E-27,  1.13091099E-26,  3.10729167E-11, -6.46234854E-27,
    -1.61558713E-26, -1.29246971E-26,  6.21458333E-11, -8.07793567E-27, -6.46234854E-27,  3.10729167E-11, -8.07793567E-27, -6.46234854E-27,  3.10729167E-11, -8.07793567E-27, -6.46234854E-27,  3.10729167E-11,
     3.10729167E-11,  1.13091099E-26, -8.07793567E-27,  6.21458333E-11,  2.26182199E-26, -1.93870456E-26,  3.10729167E-11,  1.13091099E-26, -9.69352280E-27,  3.10729167E-11,  1.13091099E-26, -9.69352280E-27,
     1.13091099E-26,  3.10729167E-11, -6.46234854E-27,  2.26182199E-26,  6.21458333E-11, -1.29246971E-26,  1.13091099E-26,  3.10729167E-11, -6.46234854E-27,  1.13091099E-26,  3.10729167E-11, -6.46234854E-27,
    -8.07793567E-27, -6.46234854E-27,  3.10729167E-11, -1.93870456E-26, -1.29246971E-26,  6.21458333E-11, -8.07793567E-27, -6.46234854E-27,  3.10729167E-11, -8.07793567E-27, -6.46234854E-27,  3.10729167E-11,
     3.10729167E-11,  1.13091099E-26, -8.07793567E-27,  3.10729167E-11,  1.13091099E-26, -8.07793567E-27,  6.21458333E-11,  2.26182199E-26, -1.93870456E-26,  3.10729167E-11,  1.13091099E-26, -9.69352280E-27,
     1.13091099E-26,  3.10729167E-11, -6.46234854E-27,  1.13091099E-26,  3.10729167E-11, -6.46234854E-27,  2.26182199E-26,  6.21458333E-11, -1.29246971E-26,  1.13091099E-26,  3.10729167E-11, -6.46234854E-27,
    -8.07793567E-27, -6.46234854E-27,  3.10729167E-11, -9.69352280E-27, -6.46234854E-27,  3.10729167E-11, -1.93870456E-26, -1.29246971E-26,  6.21458333E-11, -8.07793567E-27, -6.46234854E-27,  3.10729167E-11,
     3.10729167E-11,  1.13091099E-26, -8.07793567E-27,  3.10729167E-11,  1.13091099E-26, -8.07793567E-27,  3.10729167E-11,  1.13091099E-26, -8.07793567E-27,  6.21458333E-11,  2.26182199E-26, -1.93870456E-26,
     1.13091099E-26,  3.10729167E-11, -6.46234854E-27,  1.13091099E-26,  3.10729167E-11, -6.46234854E-27,  1.13091099E-26,  3.10729167E-11, -6.46234854E-27,  2.26182199E-26,  6.21458333E-11, -1.29246971E-26,
    -8.07793567E-27, -6.46234854E-27,  3.10729167E-11, -9.69352280E-27, -6.46234854E-27,  3.10729167E-11, -9.69352280E-27, -6.46234854E-27,  3.10729167E-11, -1.93870456E-26, -1.29246971E-26,  6.21458333E-11
    });

    // Check
    ASSERT_TRUE (coupledMTarget * factor == coupledMTest * factor);
}

TEST_F (Tetra4CalculatorTest, ComputeElementStrainStressWithCORDMeq0)
{
    // Test data
    const Matrixd TS (3, 3);
    Matrixd TMat (3, 3);
    SolidMaterialParser::ComputeMaterialCoordTrans (property.m_mCId, T.Transpose (), TS, _OUT TMat);
    SolidMaterialParser::ComputeElasticMatrix (DB_SERVICEFACTORY, property.m_mId, TMat, _OUT D, _OUT rho);

    Matrixd eleStrainTest (6, 1);
    Matrixd eleStressTest (6, 1);
    eleCalculator.ComputeElementStrainStress (localCoord, D, eleDisplacement, TMat, _OUT eleStrainTest, _OUT eleStressTest);

    // Target data
    const Matrixd eleStrainTarget (6, 1, { 1.451033E-06, -3.430588E-05,  1.126646E-04, -9.247534E-06,  2.189474E-04,  1.031818E-04});
    const Matrixd eleStressTarget (6, 1, { 9.903652E+00,  4.127536E+00,  2.786891E+01, -7.469162E-01,  1.768421E+01,  8.333916E+00});

    // Check
    ASSERT_TRUE (eleStrainTarget * factor == eleStrainTest * factor);
    ASSERT_TRUE (eleStressTarget == eleStressTest);
}

TEST_F (Tetra4CalculatorTest, ComputeElementStrainStressWithCORDMgt0)
{
    // Input data
    property.m_mCId = 100;

    // Test data
    const Matrixd TS (3, 3);
    Matrixd TMat (3, 3);
    SolidMaterialParser::ComputeMaterialCoordTrans (property.m_mCId, T.Transpose (), TS, _OUT TMat);
    SolidMaterialParser::ComputeElasticMatrix (DB_SERVICEFACTORY, property.m_mId, TMat, _OUT D, _OUT rho);

    Matrixd eleStrainTest (6, 1);
    Matrixd eleStressTest (6, 1);
    eleCalculator.ComputeElementStrainStress (localCoord, D, eleDisplacement, TMat, _OUT eleStrainTest, _OUT eleStressTest);

    // Target data
    const Matrixd eleStrainTarget (6, 1, { 3.684496E-06,  1.634322E-04, -8.730693E-05,  6.234924E-05,  1.084449E-04,  4.554705E-05});
    const Matrixd eleStressTarget (6, 1, { 1.026444E+01,  3.606983E+01, -4.434174E+00,  5.035900E+00,  8.759015E+00,  3.678800E+00});

    // Check
    ASSERT_TRUE (eleStrainTarget * factor == eleStrainTest * factor);
    ASSERT_TRUE (eleStressTarget == eleStressTest);
}

TEST_F (Tetra4CalculatorTest, ComputeElementStrainStressWithCORDMeqNegativeTwo)
{
    // Input data
    property.m_mCId = -2;

    // Test data
    Matrixd TS (3, 3);
    Tetra4ElementCalculator::ComputeInternalCoordTrans (pointSet, _OUT TS);
    Matrixd TMat (3, 3);
    SolidMaterialParser::ComputeMaterialCoordTrans (property.m_mCId, T.Transpose (), TS, _OUT TMat);
    SolidMaterialParser::ComputeElasticMatrix (DB_SERVICEFACTORY, property.m_mId, TMat, _OUT D, _OUT rho);

    Matrixd eleStrainTest (6, 1);
    Matrixd eleStressTest (6, 1);
    eleCalculator.ComputeElementStrainStress (localCoord, D, eleDisplacement, TMat, _OUT eleStrainTest, _OUT eleStressTest);

    // Target data
    const Matrixd eleStrainTarget (6, 1, { 8.594104E-05,  2.228775E-07, -6.354196E-06,  4.443859E-05,  1.168120E-04,  2.365209E-04});
    const Matrixd eleStressTarget (6, 1, { 2.355204E+01,  9.705257E+00,  8.642807E+00,  3.589271E+00,  9.434815E+00,  1.910361E+01});

    // Check
    ASSERT_TRUE (eleStrainTarget * factor == eleStrainTest * factor);
    ASSERT_TRUE (eleStressTarget == eleStressTest);
}

TEST_F (Tetra4CalculatorTest, ComputeElementStrainStressWithCORDMeqNegativeOne)
{
    // Input data
    property.m_mCId = -1;

    // Test data
    const Matrixd TS (3, 3);
    Matrixd TMat (3, 3);
    SolidMaterialParser::ComputeMaterialCoordTrans (property.m_mCId, T.Transpose (), TS, _OUT TMat);
    SolidMaterialParser::ComputeElasticMatrix (DB_SERVICEFACTORY, property.m_mId, TMat, _OUT D, _OUT rho);

    Matrixd eleStrainTest (6, 1);
    Matrixd eleStressTest (6, 1);
    eleCalculator.ComputeElementStrainStress (localCoord, D, eleDisplacement, TMat, _OUT eleStrainTest, _OUT eleStressTest);

    // Target data
    const Matrixd eleStrainTarget (6, 1, { 1.269680E-04,  3.868104E-05, -8.583937E-05,  1.767282E-04,  6.417348E-05, -3.305235E-05});
    const Matrixd eleStressTarget (6, 1, { 3.017948E+01,  1.591773E+01, -4.197106E+00,  1.427420E+01,  5.183243E+00, -2.669613E+00});

    // Check
    ASSERT_TRUE (eleStrainTarget * factor == eleStrainTest * factor);
    ASSERT_TRUE (eleStressTarget == eleStressTest);
}

TEST_F (Tetra4CalculatorTest, ComputeElementStrainStressWithCORDMeq0MAT9)
{
    // Input data
    property.m_mId = 9;

    eleDisplacement = Matrixd (12, 1, {
    0.0,                   0.0,                    0.0,
    0.0,                   0.0,                    0.0,
    0.0,                  -5.1420596569910566E-8,  0.0,
    6.194814300930412E-7,  2.742581380888817E-6,   1.50564559090072E-6
    });

    // Test data
    const Matrixd TS (3, 3);
    Matrixd TMat (3, 3);
    SolidMaterialParser::ComputeMaterialCoordTrans (property.m_mCId, T.Transpose (), TS, _OUT TMat);
    SolidMaterialParser::ComputeElasticMatrix (DB_SERVICEFACTORY, property.m_mId, TMat, _OUT D, _OUT rho);

    eleDisplacement = TEle.Transpose () * eleDisplacement;

    Matrixd eleStrainTest (6, 8);
    Matrixd eleStressTest (6, 8);
    eleCalculator.ComputeElementStrainStress (localCoord, D, eleDisplacement, TMat, _OUT eleStrainTest, _OUT eleStressTest);

    // Target data
    // Target data
    const Matrixd eleStrainTarget (6, 1, { 1.304171E-08, -5.878469E-07,  2.218846E-06, -5.205910E-08,  3.772657E-06,  9.446178E-07});
    const Matrixd eleStressTarget (6, 1, { 1.177274E+01,  4.478628E+00,  2.786349E+01,  2.412916E+00,  1.768421E+01,  8.713479E+00});

    // Check
    ASSERT_TRUE (eleStrainTarget * factor == eleStrainTest * factor);
    ASSERT_TRUE (eleStressTarget == eleStressTest);
}

TEST_F (Tetra4CalculatorTest, ComputeStressStiffness)
{
    // Input data
    localCoord = Matrixd (3, 4, {
    7.824164, 6.666667, 7.550383, 10.0,
    0.0,      0.0,      2.492317, 0.0,
    2.626982, 0.0,      0.0,      0.0
    });

    const Matrixd eleEleCoordStress (6, 1, {-30.952268454271124, -30.952268454271124, -72.221959726632619, 0.0, 0.0, 0.0});

    // Test data
    Matrixd kSigmaTest (12, 12);
    eleCalculator.ComputeStressStiffness (localCoord, eleEleCoordStress, _OUT kSigmaTest);

    // Target data
    const Matrixd kSigmaTarget (12, 12, {
    -38.0664960780,0.0000000000,0.0000000000,24.8479382529,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,13.2185578251,0.0000000000,0.0000000000,
    0.0000000000,-38.0664960780,0.0000000000,0.0000000000,24.8479382529,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,13.2185578251,0.0000000000,
    0.0000000000,0.0000000000,-38.0664960780,0.0000000000,0.0000000000,24.8479382529,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,13.2185578251,
    24.8479382529,0.0000000000,0.0000000000,-36.1405829483,0.0000000000,0.0000000000,13.3196630207,0.0000000000,0.0000000000,-2.0270183253,0.0000000000,0.0000000000,
    0.0000000000,24.8479382529,0.0000000000,0.0000000000,-36.1405829483,0.0000000000,0.0000000000,13.3196630207,0.0000000000,0.0000000000,-2.0270183253,0.0000000000,
    0.0000000000,0.0000000000,24.8479382529,0.0000000000,0.0000000000,-36.1405829483,0.0000000000,0.0000000000,13.3196630207,0.0000000000,0.0000000000,-2.0270183253,
    0.0000000000,0.0000000000,0.0000000000,13.3196630207,0.0000000000,0.0000000000,-18.1248220827,0.0000000000,0.0000000000,4.8051590620,0.0000000000,0.0000000000,
    0.0000000000,0.0000000000,0.0000000000,0.0000000000,13.3196630207,0.0000000000,0.0000000000,-18.1248220827,0.0000000000,0.0000000000,4.8051590620,0.0000000000,
    0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,13.3196630207,0.0000000000,0.0000000000,-18.1248220827,0.0000000000,0.0000000000,4.8051590620,
    13.2185578251,0.0000000000,0.0000000000,-2.0270183253,0.0000000000,0.0000000000,4.8051590620,0.0000000000,0.0000000000,-15.9966985618,0.0000000000,0.0000000000,
    0.0000000000,13.2185578251,0.0000000000,0.0000000000,-2.0270183253,0.0000000000,0.0000000000,4.8051590620,0.0000000000,0.0000000000,-15.9966985618,0.0000000000,
    0.0000000000,0.0000000000,13.2185578251,0.0000000000,0.0000000000,-2.0270183253,0.0000000000,0.0000000000,4.8051590620,0.0000000000,0.0000000000,-15.9966985618
    });

    // Check
    ASSERT_TRUE (kSigmaTarget == kSigmaTest);
}